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The self-organized criticality in Ehrenfest's historical dog-flea model is analyzed by simulating 
the underlying stochastic process. The fluctuations around the thermal equilibrium in the model 
are treated as avalanches. We show that the distributions for the fluctuation length differences at 
subsequent time steps are in the shape of a g-Gaussian (the distribution which is obtained naturally 
in the context of nonextensive statistical mechanics) if one avoids the finite size effects by increasing 
the system size. We provide a clear numerical evidence that the relation between the exponent 
t of avalanche size distribution obtained by maximum likelihood estimation and the q value of 
appropriate q-Gaussian obeys the analytical result recently introduced by Caruso et al. [Phys. Rev. 
E 75, 055f0f(R) (2007)]. This rescues the q parameter to remain as a fitting parameter and allows 
us to determine its value a priori from one of the well known exponents of such dynamical systems. 
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Introduction: The term self-organized criticality 
(SOC) was first introduced by Bak, Tang, and Wiesen- 
feld (BTW) in 1987 Q. In their well known paper, the 
so-called BTW sandpile model was used to demonstrate 
that the dynamics which gives rise to the power-law cor- 
relations seen in the non-equilibrium steady states must 
not involve any fine-tuning of parameters. Namely, sys- 
tems under their natural evolution are driven at a very 
slow rate until one of their elements reaches a thresh- 
old, i.e., statistically stationary state, and this triggers 
a burst of activity (avalanche) which occurs on a very 
short time scale. When the avalanche is over, the system 
evolves again according to the slow drive until a next 
avalanche is triggered. The activity of the system in this 
way consists of a series of avalanches. There are many 
systems where the SOC paradigm has been applied, e.g. 
earthquakes, noise with 1// power spectrum, brain ac- 
tivity, river networks, biological evolution of interacting 
species, traffic jams etc. 0). 

Following the BTW sandpile model a great variety of 
models from the deterministic and stochastic to the dis- 
sipative and conservative have been introduced which ex- 
hibit the phenomenon of SOC (for an overview, see Q 
and references therein). In 1996, a random neighbor ver- 
sion of the original BTW sandpile model was presented 
by Flyvbjerg In this work, it was emphasized that a 
self-organized critical system is a driven, dissipative sys- 
tem consisting of a medium (sandpile) which has distur- 
bance propagating through it, causing a modification of 
the medium, such that eventually the medium is in a crit- 
ical state, and the medium is modified no more. More- 
over, it was shown by way of random neighbor sandpile 
model that a dynamical system with only two degrees of 
freedom can be self-organized critical and as it is the case 
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in fluctuation phenomena, the dynamics is described by 
a master equation which can be partially solved analyti- 
cally. 

Soon after Flyvbjerg's work Nagler et al. studied the 
conservative variant of random neighbor sandpile model 
which is neither extended nor dissipative with regard to 
the amount of sand in the system but still shows SOC 
with nontrivial exponents [H, Q. This kind of analy- 
sis is not restricted to nonspatial systems and available 
also for spatial systems like one-dimensional cellular au- 
tomata [7j. The dynamics of the model described by 
Nagler et al. is given on a Fokker-Planck equation by 
introducing appropriate scaling variables. The avalanche 
size distribution which is readily obtained by solving the 
Fokker-Planck equation at an absorbing boundary ex- 
hibits a power-law regime followed by an exponential 
tail. Their model is an adaptation of the famous dog-flea 
model introduced by Ehrenfest in 1907 @. This model 
can be considered as a zero-dimensional nonspatial proto- 
type SOC model and its dynamics is different from most 
of the standard SOC models which are iV-dimensional 
spatial systems. 

The dog-flea model is a simple but typical example 
of generation-recombination Markov chain Q describing 
the process of approaching an equilibrium state in a large 
set of uncoupled two state systems together with fluctua- 
tions (avalanches) around this state. For an even number 
of states, the transition probability of fluctuations of the 
discrete time version was calculated by Kac [l(| (see also 
fill]). An identification of the model as a random walk on 
a Bethe lattice is studied in Ref. [HI]. Furthermore, it has 
recently been shown that the dog-flea model, formulated 
as a continuous time Markov chain, is a representation of 
a spin in a magnetic field fl3| . Such a representation is 
used to estimate the blocking temperature in molecular 
nano-magnets (l4| . 

In this work, we will be analyzing the SOC in the dog- 
flea model through simulation of the underlying stochas- 
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tic process that describes the natural evolution of the 
model. The analysis method that we use has recently 
been presented by Caruso et al. to interpret the SOC 
in the limited number of earthquakes (up to 689 000) 
taken from World and Northern California catalogs for 
the periods 2001-2006 and 1966-2006, respectively 
Using the same line of thought, it is our aim to ana- 
lyze the SOC feature of the dog-flea model through the 
time series of the fluctuation length. The simplicity of 
the dynamics of the dog-flea model enables us to obtain 
a large number of fluctuations for different system sizes 
in a reasonable computing time (i.e., we consider up to 
2 x 10 9 fluctuations). Thus, the obtained critical expo- 
nents for the model are very precise as it will be discussed 
in coming sections. This analysis enables us to accom- 
plish our main task, which is to provide the first rigorous 
numerical example where the relationship, proposed by 
Caruso et al., between the exponent r of avalanche size 
distribution and the q value of appropriate g-Gaussian 
(the distribution which is obtained naturally in the con- 
text of nonextensive statistical mechanics) [16[ ■ This will 
be very appealing also from nonextensive statistical me- 
chanics point of view since this treatment makes the q 
parameter to be determined a priori, which is a situa- 
tion achieved rarely up to now. 

The model and numerical procedure: The dynamics of 
the dog-flea model has simple rules. The model has iV 
dynamical sites represented by the total number of fleas 
shared by two dogs (dog A and dog B). Suppose that 
there are Na fleas on dog A and Nb fleas on dog B 
leading to a population of fleas N = Na + Nb- For 
convenience, N is assumed to be even. In every time 
step, a randomly chosen flea jumps from one dog to the 
other. Thus, we have N A -> N A ± 1 and N B -> N B =F 1- 
The procedure is repeated for an arbitrary number of 
times. In long time run, the mean number of fleas on 
both dog A and dog B converges to the equilibrium value, 
(Na) = (N B ) = N/2 with the fluctuations around it. A 
single fluctuation is described as a process that starts 
once the number of fleas on one of the dogs becomes 
larger (or smaller) than the equilibrium value N/2 and 
stops when it gets back to it for the first time. Thus, the 
end of one fluctuation specifies the start of the subsequent 
one. The length (A) of a fluctuation is determined by the 
number of time steps elapsed until the fluctuation ends. 

It is straightforward to obtain the master equation of 
the process that describes the time evolution of the prob- 
ability to find a specified number of fleas on one of the 
dogs. Assuming that after t steps there are iV^t) = I 
fleas on dog A, at the subsequent time step there are 
only two possibilities, f^l + lorl^l-1 with 
the transition probabilities W(£ + 1\£) = (N - l)/N and 
W(£ — l\£) — £/N, respectively. Then, the time evolution 
of the probability P(£, t) to find £ fleas on dog A at time 
t obeys the following master equation, 

£+1 N — £+ 1 
P(£, t+1) = ~i^P(t+l, t)+ -— P(£-l, t). (1) 



Introducing appropriate scaling variables Eq. ([I]) can be 
written in the form of a Fokker-Planck equation by which 
the fluctuation distribution is reviewed analytically 

Distribution of fluctuation length and returns: As 
it was first demonstrated by BTW sandpile model, a 
generic signature of SOC is the presence of a power- 
law as well as finite size scaling in the size or the dura- 
tion distribution of the avalanches. Recently, a power- 
law regime following an exponential tail in the fluc- 
tuation length distribution for the Ehrenfest's dog-flea 
model has been reported for a very limited system size 
(i.e., N — 2500) @. In our paper, in order to ana- 
lyze the SOC in the dog-flea model through the fluctu- 
ation length distribution we simulate the corresponding 
stochastic process for seven different values of N namely, 
N = 10 2 , 10 3 , 5 x 10 3 , 10 4 , 10 5 , 10 6 , and 10 7 . For con- 
venience, let us group the first four different system sizes 
as " small Ns" and the remaining sizes as "large TVs". In 
Fig. Q^a) and (b) we plot the distribution of the fluctua- 
tion length time-series X(t) for the small TVs and large TVs, 
respectively. In order to have good statistics 10 9 fluctua- 
tions for the small iVs group and 2 x 10 9 fluctuations for 
the large iVs group have been considered. In both cases 
the fluctuation distributions have a power-law regime, 
P(X) ~ A~ T while in the small Ns group the power-law 
regime is followed by an exponential decay because of 
the finite-size effect. For the small As group one can 
control if the fluctuation length distribution P(X) obeys 
the following finite size scaling behavior, 

P ^ ~ J^f (jk) • ( 2 ) 

where / is a suitable scaling function and 7 and £ are 
critical exponents describing the scaling of the distribu- 
tion function. In the inset of Fig. [Ha), a clear data 
collapse of P(X) is shown for the small Ns group (i.e., 
N = 10 2 , 10 3 , 5 x 10 3 , and 10 4 ). This data collapse in- 
dicates that the fluctuation length distributions of small 
As satisfy the finite size scaling hypothesis very well. 
The obtained critical exponents are 7 ~ 1.517 and £ = 1. 
As it is seen from Fig. [IJb), these values of critical ex- 
ponents are in agreement with the finite size scaling hy- 
pothesis since for asymptotically large N, P(X) ~ A~ T 
with t — j/( ~ 1.517. The value of r is obtained by the 
maximum likelihood estimation (MLE) and this method 
enables us to determine this exponent of the model as 
accurate as ±1.156 x 10" 5 (l7j . 

Now we are at the position to introduce the dis- 
tribution of returns, i.e., the differences between fluc- 
tuation lengths obtained at consecutive time steps, as 
AX(t) = X(t + 1) — X(t). It should also be noted that, in 
order to have zero mean, the returns are normalized by 
introducing the variable x as 

a: = AA-(AA), (3) 

where (■ ■ • ) stays for the mean value of the given data 
set. The signal of the distribution of returns reveals very 
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FIG. 1: (color online) Fluctuation length distributions for the small iVs and for the large JVs groups are given in (a) and (b), 
respectively. In the inset of (a), we also present data collapse of finite size scaling given in Eq. ((2)1 for small iVs group. The 
critical exponents derived from the fit are 7 ~ 1.517 and f = 1. The full black line in (b) represents the fitting curve of the 
distribution with slope r ~ 1.517 which has been obtained by maximum likelihood estimation. The distributions have an 
arbitrary normalization such that P(A = !) = !. 
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FIG. 2: (color online) The distributions of returns, i.e., the fluctuation length differences AA(t) = A(t + 1) — A(t), normalized 
by introducing the variable x = AA — (AA) are shown in (a) for the small Ns group and in (b) for the large iVs group. For 
comparison, standard Gaussian and q-Gaussian curves are drawn by black dashed and full lines, respectively. See text for 
further details. In insets, the central parts of the distributions are emphasized. 



interesting results on the criticality of the dog-flea model. 
This approach is used in recent studies on turbulence fhSj 
and the time-series of real earthquakes [l5l |. 

In Fig. [21 we plot the distribution of the returns AA(t) 
obtained from 10 9 fluctuations for each different system 
sizes in the small Ns group (a), whereas in the group of 
large Ns (b) 2 x 10 9 fluctuations are considered. What is 
common for both cases is that none of them has return 
distributions which can be approached by a Gaussian. As 
the system size N increases, leading to a longer power-law 
regime in the fluctuation length distribution, the return 
distribution curves become to exhibit a convergence to 
a kind of fat tailed distribution. When the system size 
is large enough, the exponential decay of the fluctuation 



length distribution (see Fig. QJb)) is postponed to larger 
sizes and the finite size effects get invisible up to more 
than four decades. In this case the distribution of the 
returns can be fitted by a g-Gaussian given by 

P(x) = P(0)[(l + P(q-l)x 2 } 1 ^\ (4) 

where (3 characterizes the width of the distribution and 
q is the index of nonextensive statistical mechanics [Tf3 | 
(black full lines in Figs.^a) and (b)). In Eq. Q, q ^ 1 
indicates a departure from the Gaussian shape while nor- 
mal Gaussian distribution can be recovered again in the 
q — > 1 limit. Here, it is worth mentioning that our re- 
sults in Fig. [5] clearly show the connection between crit- 
icality and the appearance of g-Gaussian, namely, wider 
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FIG. 3: (color online) (a) Distribution of returns for a repre- 
sentative case of large As group (N — 10 6 ) is given by full 
green circles. The (/-Gaussian curve with q — 2.35 and /3 = 35 
is shown by full black line. This value of the q is obtained 
by substituting r = 1.517 into Eq. iJBJ. A standard Gaussian 
curve is drawn by dashed black line for comparison. In the 
inset, the central part of the distribution is given in order to 
emphasize that the distribution approaches almost perfectly 
to the (/-Gaussian not only in the tails but also in the center, 
(b) In order to better visualize how well the used (/-Gaussian 
approaches to the distribution, we plot the same P(x) versus 
1 + P(q — l)x 2 . A straight line with a slope 1/(1 — q) is ex- 
pected for a perfectly (/-Gaussian shaped distribution. Data 
points (green circles) and the slope with q — 2.35 (black line) 
constitute a clear evidence towards this tendency. 



the critical regime persists, longer the tails of returns 
distribution follow {/-Gaussian. This kind of interpreta- 
tion might also be useful in understanding the difference 
between two recent experimental works on velocity dis- 
tributions in optical lattices [lij In pjj], velocity 
distributions are found to approach a double-Gaussian 
shape, whereas in (20| they are reported to converge to 
a g-Gaussian. The reason for this discrepancy seen in 
the results of essentially the same experiment might be 
that in the latter the system may be set exactly at the 
criticality, whereas in the former it is not. 

At this point, we should recall the important result 



reported by Caruso et al. [15| relating the r exponent of 
the avalanche size distribution with the q parameter of 
the g-Gaussian. As it was emphasized in their work, if 
there is no correlation between the size of two events, the 
probability of obtaining the difference AA = X(t+6) — X(t) 
(5 is an integer describing the correlation length and in 
our case 5 = 1) is given by 



P(AA) = K- 



"(2T-1) 

2r - 1 



F x t,2t-1;2t;-J 



IAAI 



(5) 



where K is a normalization factor, e is a small posi- 
tive value and 2-F1 is the hypergeometric function. The 
curve of this r dependent probability density function 
P(AA) can be approached by means of g-Gaussian with 
e- independent q value. In Ref. [HI], by evaluating Eq. §5§ 
for various values of r, a relation between the power-law 
exponent r and q is reported as 



,1.19t" 



(6) 



Although this relation is obtained in [l5| by Caruso et 
al., they could not check its validity since the earthquake 
data that they analyzed was not adequate to obtain the 
r value with high precision. Consequently, they still used 
q parameter as a fitting parameter. On the other hand, 
since the power-law exponent is very accurate in our case, 
we can substitute its value (r = 1.517) obtained by MLE 
into Eq. ([6]) which gives the q value as q = 2.35. This 
value is obviously the one that we should use in the q- 
Gaussian to check whether the return distribution can be 
approached by this. It is worth mentioning here that the 
q parameter is not a fitting parameter anymore. In Fig. [2] 
we also include this result together with a Gaussian curve 
for comparison. It is clear that, for very small Ns, the 
convergence to g-Gaussian is only in the central part (see 
the inset of Fig.[2ja)), whereas it develops more and more 
towards the tails as N increases. Eventually, for large 
enough A^s for which finite size effects are invisible inside 
the obtained region, the g-Gaussian curve is perfectly 
approached including the center and tails. 

In order to further strengthen our results, we consider 
one of the appropriate system size (N = 10 6 ) separately 
in Fig. [3l A very clear convergence of the return distri- 
bution to the g-Gaussian can be seen everywhere for the 
available data (including the very central part, see the 
inset of Fig. [21a)). Moreover, to check how well the ob- 
tained g-Gaussian curve approaches the returns distribu- 
tion, a log-log plot of Eq. (fj| is given in Fig.[3]Jb). A per- 
fect straight line with the slope 1/(1 — g) is the expected 
behavior for this type of representation if the curve is an 
exact g-Gaussian and as it is seen very clearly, the be- 
havior of the return distribution fulfills this tendency ex- 
hibiting a seven decade power-law with the slope 1 / ( 1 — g) 
which gives the already obtained q value, q — 2.35. 

Conclusion: We analyze the SOC in the Ehrenfest's 
dog-flea model through the probability distributions 
of the fluctuation length (avalanche size distributions) 
and of the differences between the fluctuation lengths 



5 



at subsequent time steps (returns distributions) by 
simulating the stochastic process of the model. Our 
extensive simulations enable us to determine the power- 
law exponent r of the avalanche size distribution with 
an extreme precision. Then, the behavior of the returns 
distributions is analyzed and numerically shown that it 
converges to a g-Gaussian with q = 2.35, a value coming 
directly (and a priori) from Eq. ^ which makes q 
parameter to be related to one of the well known power- 
law exponents of such model systems (which means 
that q is not a fitting parameter anymore). This is the 
main result of the present letter and important from 
(at least) three point of view: (i) this constitutes the 
first reliable verification of Caruso et al. relation since, 
due to insufficient data set of earthquakes, they were 
unable to provide a clear evidence for their own relation; 
(ii) this result is achieved using a simple, prototype 
SOC model (different from the one used by Caruso et 
al.) which can be considered as the first clue on the 
generality of these results rather than being specific only 
to this model; (iii) this treatment makes the q parameter 
of the g-Gaussian to be determined a priori which 
constitutes a rather rare achievement in the literature 
due to technical difficulties. ^From the analysis of return 



distributions from small TVs to large iVs, it is shown that 
the convergence to appropriate g-Gaussian starts from 
the central part and gradually develops towards the tails 
as N increases. This is a kind of expected behavior 
since, from our simulations it is also evident that the 
power-law regimes of the avalanche size distributions 
for small Ns are followed by exponential decays due to 
finite size effects and this obviously deteriorates the true 
behavior. Of course, for large enough iVs, this effect is 
postponed further and further to avalanche sizes that are 
not inside the region we are considering. Moreover, one 
could conclude that, as N — ► oo the power-law regime 
of avalanche size distribution is expected to continue 
forever, then the corresponding return distribution 
appears to converge to the g-Gaussian for the entire 
region. Finally, it is worth to mention that the behavior 
observed and reported here for the zero dimensional 
prototype SOC model of Ehrcnfest is by no means spe- 
cific and limited to this model, but seems to appear as a 
rather common phenomenon for several SOC models [2l| . 
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